An Analysis of Polynomial Composition Algorithms

An analysis is made of the runtime of a previously published algorithm for polynomial composition Two new more e cient algorithms are presented One of these algorithms is optimal while the other algorithm is numerically more stable than the optimal one Additionally as a generalization of polynomial composition we show how to compose a multia ne function with a set of polynomials as an extension to an earlier algorithm for com posing two polynomial functions With this extension we are able to perform degree raising with composition Introduction Many CAGD operations can be implemented using polynomial composition As discussed in we can implement evaluation subdivision freeform deformation conversion between tensor product and B ezier simplex forms using polynomial composition In that paper an e cient algorithm for composing two polynomials in B ezier simplex form was presented No analysis was given for the runtime of this algorithm however In Section we analyze the runtime of this composition algorithm In Section we then present a variation on this algorithm that is more e cient A second variation that is optimal is presented in Section However this optimal algorithm is numerically less stable than the non optimal one motivating the presentation of both techniques Further while we can recast many CAGD problems as polynomial composition we have been unable to implement degree raising using polynomial composition That is given the control points for a degree n polynomial function F nd the control points for the degree n representation of this function In Section we present a generalization of the composition algorithm for two B ezier simplices and show how it can be used to degree raise a B ezier simplex All the pseudo code used by the algorithms appears in Appendix B Notation and Background In this section we brie y summarize the basic concepts of blossoming and we introduce a number of notational conventions that have been invented to simplify the manipulations needed in Section University of Waterloo CS The notation used in this paper is a simple extension of the notation in the DeRose et al paper and is summarized in Table Multi indices and Hyperindices Most of our notation deals with the indexing of multivariate polynomials While univariate Bern stein polynomials are commonly indexed with integers multivariate Bernstein polynomials are most easily indexed using tuples of non negative integers such as i ik i ik We shall refer to such tuples as multi indices We denote the norm of a multi index by j j and de ne it to be the sum of the components of We use the symbol IIk to stand for the set of all multi indices i ik with j j d For example II f g The symbol ej will be used to denote a multi index that is in the jth component and zero in the rest of the components denotes a multi index whose components are all zero Addition and subtraction of multi indices is de ned componentwise We shall also have occasion to use tuples of multi indices such as I m II d m k where II d m k denotes the m fold Cartesian product II d k II dm k for notational convince we ignore d We shall refer to such tuples as hyper indices When d dm d we denote II d m k as II d m k Again we use the notation jI j to denote the sum of the components of I Notice that whereas j j is an integer jI j is a multi index Therefore for I II d m k the expression jjI jj evaluates to the integer j j j mj Note that for most of the paper we will use the II d m k notation We will not use the II d m k notation until Section Bernsteins and B eziers Using this indexing notation the k variate Bernstein polynomials of degree d can be de ned by B b bk d b b i b ik k where i ik II d k b bk IR and where d d i i ik is the multinomial coe cient Notice that each of the Bernstein polynomials is a homogeneous polynomial of degree d It is known that for every polynomial Q X Y of degree d where X is an a ne space of dimension k and Y is an a ne space of arbitrary dimension there exist unique points fV g II k in Y such that Q u X II k V B d b bk An Analysis of Polynomial Composition Algorithms where b bk are the barycentric coordinates of u X relative to a simplex x xk of points in X For a description of a ne spaces simplices and barycentric coordinates see for instance Farin or DeRose A polynomial Q when expressed as in Equation is called a B ezier simplex the points V are called the control net of Q with respect to the domain simplex Note that a degree d B ezier simplex with a domain of dimension k has jIId j d k k control points We shall often writeB u with the understanding that u should be replaced with its barycentric coordinates relative to the appropriate domain simplex If I m II d m k we de ne B d m I u to be the product B u Bm m u We shall make use of the following product relation satis ed by the Bernstein polynomials B d m I u C I B j dj jIj u where C I is a combinatorial constant given by C I j j j mj m jjI jj jI j This relation is easily proved using simple manipulation of the Bernstein polynomials c f Blossoms Ramshaw discovered how to exploit a connection between B ezier simplices and symmetric multi a ne maps A map q u ud is said to be multi a ne if it is a ne when all but one of its arguments are held xed it is said to be symmetric if its value does not depend on the ordering of the arguments Associated with every polynomial Q X Y of degree d there is a unique symmetric d a ne map q X d Y that agrees with Q on its diagonal the diagonal of a multi a ne map q u ud is the function obtained when all arguments are equal Q u q u u u Ramshaw refers to this multi a ne map q as the blossom of Q Ramshaw has shown that the B ezier control net for a polynomial relative to a domain simplex can be obtained by evaluating the polynomial s blossom at the vertices of the simplex More precisely if X is an a ne space of dimension k Q X Y a polynomial of degree d whose blossom is q and x xk a simplex in X then Ramshaw has shown that the B ezier control net of Q relative to is given by V q d z x x z i x x z i xk xk z ik for all i ik II d k Using Farin s blossom notation this equation can be written as V q x i x i x ik k University of Waterloo CS Polynomial Composition The basic polynomial composition problem is the following Given A ne spaces X Y and Z of dimensions KX KY and KZ respectively control points fG g II KX de ning a B ezier simplex G X Y of degree relative to a domain simplex X X and control points fF pg p IImKY de ning a B ezier simplex F Y Z relative to a domain simplex Y Y Find The control points fH g II KX of the degreem B ezier simplexH F G relative to X Solution If f denotes the blossom of F then H X I II m KX jI j C I f GI II L KX where GI with I m is an abbreviation for G G m A proof of this result can be found in the paper by DeRose et al and a generalization of this proof can be found in Section In the following sections we will make frequent use of the number of G s control points which we will denote as G Note that G KX KX with G KX and de ned above Implementation and Analysis DeRose et al used symmetry and reuse of intermediate results to reduce the computation cost of computing the control points ofH using Equation In this section we discuss these two techniques and analyze their runtime behavior Note that the code in Figures and appeared in Further note that unless otherwise stated our algorithms will evaluate the blossom f GI one argument at a time in the order given by I de Casteljau s Algorithm B ezier curves and surfaces are commonly evaluated using de Casteljau s algorithm This algorithm uses repeated linear interpolation to evaluate B ezier simplices In addition to evaluating the B ezier simplex for position the intermediate results can be used to compute derivatives Further the algorithm can be extended to evaluate a single argument of the blossom See Farin s book for more details on de Casteljau s algorithm Note that the intermediate values used in a de Casteljau evaluation of a degree d dimension k B ezier simplex can be stored in simplex of one higher dimension i e with d k k storage Further each point in this simplex other than the initial control points is computed using one linear combination Thus a full de Casteljau evaluation requires d k k linear combinations during its computation An Analysis of Polynomial Composition Algorithms Notation Description IR The set of real numbers IR n dimensional real space a z Integers and real numbers X Y etc A ne spaces KX The dimension of the a ne space X a z A ne points Domain simplex a z Multi indices j j i ik where i ik IIk The set of all multi indices i ik with j j d I J I etc Hyper indices II m k The set of all hyper indices I m with m II d k II d m k The set of all hyper indices I m with j II dj k jI j m where I m jjI jj dm where I II m k j j The multinomial coe cient j j i ik i ik B u The th Bernstein polynomial of degree d B d m I B d u B m u I m C I A combinatorial constant given by j j j m j m jjIjj jIj q u ud Blossom of a B ezier simplex Q u of degree d f GI f G G d with I II d m k G KX KX the number of control points of G Table Summary of notation University of Waterloo CS In the initial analysis of runtime we will count the number of linear combinations performed by the composition algorithms A more complete runtime analysis can be found in Section where we show that the cost of the more sophisticated algorithms is proportional to the number of linear combinations performed In Figure we present some code for performing partial blossom evaluations All the linear combinations are performed in EvalBlossomArgument which evaluates one blossom argument by performing a single step of de Casteljau s algorithm Naive Algorithm A naive implementation of the composition algorithm would compute each H by performing a full evaluation of f for every possible hyperindex I II m KX where jI j Note that this means we will be evaluating f exactly once for each hyperindex I As mentioned earlier a full de Casteljau evaluation requires m KY KY linear combinations Since the number of multi indices of dimension KX and degree is KX KX G the number of hyper indices will be G Thus this algorithm requires m KY KY G linear combinations We shall refer to this algorithm as the Naive Algorithm One Permutation Algorithm The rst way to improve the Naive Algorithm is to use the symmetry property of blossoms This property allows us to evaluate f only at hyperindices that are unique up to permutations Using such hyperindices we only care how many times each multi index j appears in I not where each j appears Thus we can convert the hyperindex I to a multi index where ij is the number of times j appears in I The dimension of is one less than the number ofG s control points and the degree is the number of multi indices in the hyperindex I which in the case we are considering is m Thus II G For example if G has ve control points G G and if F is fourth degree then the hyperindex I can be thought of as the multi index II Normally we prefer to think of hyperindices and multi indices as di erent since they are used in di erent ways and if we care about the order of the multi indices within the hyperindex then the above correspondence fails to hold However such a correspondence is useful in counting the number of hyperindices that are unique up to permutations which is